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Abstract 

We present a new approximation to the normal distribution quan- 
tile function. It has a similar form to the approximation of Beasley 
and Springer [3J, providing a maximum absolute error of less than 
2.5 • 10~ 5 . This is less accurate than [3], but still sufficient for many 
applications. However it is faster than [3]. This is its primary bene- 
fit, which can be crucial to many applications, including in financial 
markets. 



1 Introduction 

The use of the inverse of the CDF for a probability distribution, also known as 
the quantile function, is widespread in statistical modelling (see, for example, 

BED- 

During recent work, the need arose for a fast and reasonably accurate ap- 
proximation to the normal distribution quantile function, N" 1 ^). Accuracy 
similar to the approximation in Equation 26.2.23 of pQ was sufficient (max 
absolute error less than 4.5 • 10~ 4 ). But speed was crucial. 

The approximation of Beasley and Springer [3], along with related ap- 
proximations such as Acklam's [2] , provides improvements in terms of both 
accuracy and speed. 

Both the Acklam and the Beasley-Springer approximations are based on 
the same ideas: 

(1) consider narrow tails separately from a wide central area 

(2) use a rational function of x to approximate N~ 1 (x) in this wide central 
area (avoiding expensive operations like log and sqrt) 

(3) take advantage of the fact that N~ 1 (x — 1/2) is an odd function. 
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The second and third ideas suggest that for the central region, we consider 
rational approximations of the form 



where F is a rational function. The approximations of Acklam, Beasley- 
Springer, and others for the central region are of this form. 

The Beasley- Springer approximation for the central region is sometimes 
called a (3,4) scheme, since the numerator of F is cubic in (x — 1/2) 2 and 
the denominator of F is of degree 4 in (x — 1/2) 2 . Similarly, the Acklam 
approximation is called a (5, 5) scheme. 

2 New Approximations 

For increased speed, here we consider a (2,2) scheme for the central region 
and a (3, 2) scheme for the tails. 

We chose the boundaries between the central region and the tails to be 
at 0.0465 and 0.9535, since with the above schemes and boundaries the max- 
imum absolute error in both regions was nearly the same and both slightly 
less than 2.5 • 10~ 5 . 

2.1 Central Region 
2.1.1 0.0465 < p < 0.9535 

Put q = p - 0.5 and let r = q 2 . For 0.0465 < p < 0.9535, define 



{x-l/2)F{{x-l/2)\ 




where 



0.389422403767615, 

-1.699385796345221, 

1.246899760652504, 

0.195740115269792, 

-0.652871358365296, 

0.155331081623168, 



-0.839293158122257. 
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The benefit of the second expression is that we save one multiplication 
by using it. Similarly, normalising the denominator so that the leading co- 
efficient is 1, rather than the constant coefficient as some authors do, also 
saves another multiplication. 

There are 12 points of maximum error (also known as alternating points) 
in the interval [0.0465,0.9535]: 



[jp, err abs ) 


(p, err abs ) 


(0.046500,2.494327- 10" 5 ) 


(0.592289,2.494326- HT 5 ) 


(0.054264,2.494331 • 10" 5 ) 


(0.752182,2.494327- 10" 5 ) 


(0.081621,2.494328- HT 5 ) 


(0.859308,2.494323- 10" 5 ) 


(0.140694,2.494323- 10~ 5 ) 


(0.918381,2.494328- HT 5 ) 


(0.247820,2.494327- 10~ 5 ) 


(0.945738,2.494331 • 10" 5 ) 


(0.407712,2.494326- HT 5 ) 


(0.945350,2.494327- HT 5 ) 



From the theorems of Chebyshev and de la Vallee Poussin (see [4J Sec- 
tion 5.5]), it follows that f ce ntrai{p) is essentially the best possible rational 
approximation of (2, 2) scheme. 

For comparison, the maximum absolute error of the "central" approxi- 
mation in [3] is under 1.85 • 10~ 9 . 

This approximation was found using the minimax function within the 
numapprox package of Maple: 

Digits: =60: with (numapprox) : 
uBnd: =0.4535" 2: 

minimax (x->inverseCDFCentralRatApprox(x) ,0. .uBnd, [2,2] ,x->sqrt(x)) ; 

where 

inverseCDFCentralRatApprox(x) is the function A^~ 1 (y / x + 1/2) /y/x, 
uBnd is the range we want the approximation over, 

[2, 2] specifies that we want the degree of both the numerator and the de- 
nominator to be 2, and 

\fx is the weight function we use, since we want to get the best approximation 
to iV~ 1 ( v /i+ 1/2) rather than N~ 1 ( y /x + l/2)/sfx. 

We tried other values of uBnd near 0.4535, but the smallest maximum 
absolute error was found with this particular value. 
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2.1.2 0.025 < p < 0.975 

The use of an even wider central region may be preferred, as this can provide 
further performance gains by reducing the expensive log and sqrt operations 
required for the tails. 

We give one such example here (found as above using Maple, but with 
uBnd=0.475). 

Put q = p - 0.5 and let r = q 2 . For 0.025 < p < 0.975, define 



a = 0.151015505647689, 

Ql = -.5303572634357367, 

a 2 = 1.365020122861334, 

b = 0.132089632343748, 

^ = -.7607324991323768. 



The maximum absolute error for this approximation is less than 1.16- 10~ 4 
which occurs near p = 0.9692. While this error is much larger than the error 
in the previous section, it is still well smaller than the maximum error for 
the Abramowitz-Stegun approximation (4.5 • 10~ 4 ). 




where 



2.2 Tails 



2.2.1 e" 372 / 2 < p < 0.0465 



For 5.3 ... • 10~ 298 = e- 3j2/2 <p< 0.0465, put r = ^\og(l/p 2 ) and define 




where 



c = 16.896201479841517652, 
ci = -2.793522347562718412, 
c 2 = -8.731478129786263127, 
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c 3 = -1.000182518730158122, 

c = 16.682320830719986527, 

c[ = 4.120411523939115059, 

4 = 0.029814187308200211, 

d = 7.173787663925508066, 

d l = 8.759693508958633869. 

As with the "central" approximation, this approximation was also found 
using the minimax function within the numapprox package of Maple: 

Digits: =60: with (numapprox) : 

v:=0.0465: 

uBnd: =0.4535"2: 

minimax (y->inverseCDF (exp (-y*y/2) ) , sqrt (log(l/v~2) ) . .37, [3,2] ) ; 

Note that since we are approximating N~ 1 (x) itself here, we do not include a 
weight function in the arguments of the minimax function and so the default 
weight function 1 is used. 

The maximum absolute error in this case is less than 2.458 ■ 10 -5 . 

2.2.2 0.9535 < p < 1 - e" 37 '/ 2 

Due to the symmetry of N~ 1 (p) about p — 1/2, we approximate iV _1 (p) by 
-ftauO- - p) (note that here r = 0og(l/(l - p) 2 )). 



3 Abramowitz and Stegun Approximations 

Having found the above new approximations, we turned our attention to the 
approximations in Equations 26.2.22 and 26.2.23 of [I]. As those authors 
note, these approximations are from [5]. In particular, Sheets 67 and 68 on 
pages 191-192 of [BJ. 

If we restrict our attention to ranges like e~ 3j2 / 2 < p < 1 — e~ 372//2 (this 
includes almost the entire IEEE- 754 range of representable real numbers), 
then we can improve on the approximations of Abramowitz and Stegun. 

For example, in this range, we can replace Equation 26.2.23 of pQ with 

_ c 2 t 2 + Cjt + c 

Xp ~ d^ + d 2 t 2 + d l t + l + e[ph 
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where \e(p)\ < 8 • 10 5 and 



d 2 
d 3 



d 



c = 



Cl = 



c 2 = 



1 



2.653962002601684482, 
1.561533700212080345, 
0.061146735765196993, 
1.904875182836498708, 
0.454055536444233510, 
0.009547745327068945. 



This is over five times more accurate than the approximation in [T] . How- 
ever, as one increases the range even closer to and 1, the max absolute 
increases until we obtain Equation 26.2.23 of pp. The near-best possible na- 
ture of Equation 26.2.23 is illustrated by the graph in Sheet 68 of [6j showing 
that Chebyshev's theorem nearly holds for this approximation. 

Note also that this approximation shows the justification for the use of 



log(l/p 2 ) in these tail approximations. As p — > 0, iV 1 (p) approaches 



4 Performance 

Using Java (JDK 1.6.0J.7), we coded the following approximations in order 
to compare their performance. 

• the Abramowitz-Stegun approximation (AS in the table below) 

• the Beasley- Springer approximation (BS in the table below) 

• the approximation from Section 2 using the central region approximation 
in Section 2.1.1 (Rat22A in the table below) 

• the approximation from Section 2 using the central region approximation 
in Section 2.1.2 (Rat22B in the table below). 

In each case, we calculated the approximation 200,000 times for each p 
from 0.001 to 0.999 with 0.001 as our step size. These calculations were done 
on a Dell Inspiron 1525, running Windows Vista and using an Intel Core 2 
Duo T5800 2.00 GHz CPU. The times in milliseconds for each approximation 
are given in the table below. 





log(l/p 2 ) plus a quantity that approaches as p does. 
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method 


time(ms) 


AS 


25,210 


BS 


10,212 


Rat22A 


8052 


Rat22B 


6649 



As one would expect, the new approximations given here are faster than 
the currently known ones. The comparison between Rat22A and Rat22B is 
also interesting, as it shows the impact of the calculation of the log and sqrt 
operations. Although these operations only need to be performed for a small 
subset of all values of p, reducing the number of these operations by just 
under 50% reduced the CPU time required by nearly 20%. 
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